### This script will intersect ASCII files with XY data

### Load library

library('SDMTools')

# Set in, out, and work directory

in.dir='/home1/99/jc152199/Microclimate Statistical Downscale/To Analyse/'
out.dir='/home1/99/jc152199/Microclimate Statistical Downscale/To Analyse/'
setwd(in.dir)

# Read in climate data

c.data = read.csv('solar_regress_all_days.csv',header=T)
c.data$coastdist = NA
c.data$EVI = NA

# Read in loc data

in.dir = '/home1/99/jc152199/Microclimate Statistical Downscale/Location Data/'
setwd(in.dir)
locs = read.csv('Locs.csv', header=T)
locs = locs[c(1:32,34:37),] # Remove lines 33 & 38 from 'locs' because these will contain NA for the topographic data
pnts=data.frame(east=locs$east, north=locs$north)

# Define ASCII directory

ASCII.dir = '/home1/99/jc152199/Microclimate Statistical Downscale/Topographic Data/ASCII/'
ASCII.list = c("coastdist.asc","EVI.asc")

# Read in ASCII files

coastdist.asc = read.asc(paste(ASCII.dir,ASCII.list[1],sep=""))
EVI.asc = read.asc(paste(ASCII.dir,ASCII.list[2],sep=""))

# Extract data

coastdist.data = extract.data(pnts, coastdist.asc)
EVI.data = extract.data(pnts, EVI.asc)

# Row bind locs to extracted data as a single data frame

asc.data = data.frame(site=locs$site, coastdist=coastdist.data, EVI=EVI.data)

# Round data

asc.data$coastdist = round(asc.data$coastdist,0)
asc.data$EVI = round(asc.data$EVI,4)

# Remove site '32004' from c.data, produces NA's when populating c.data$EVI

c.data = c.data[-which(c.data$site=='32004'),]

# Merge topographic data back with climate data

for (ii in 1:nrow(c.data)) {
  cat(ii,'...')             
  row.in.asc = which(asc.data$site==c.data$site[ii]) # This command identifies the row of the solar data that we're interested in
  c.data$EVI[ii] = asc.data[row.in.asc,3] # This command replaces the 'NA' value from c.data$EVI with the EVI value from EVI.data
  c.data$coastdist[ii] = asc.data[row.in.asc,2] # This command replaces the 'NA' value from c.data$coastdist with the coastdist value coastdist.data
  
}

# Write out c.data and perform linear regressions on this data set

write.csv(x=c.data, file=paste(out.dir,"microMACROtopoVEG.csv",sep=""),row.names=F)